Mayo 17, 2017

Problema módelo

Se mide una variable en una población que viene de dos subpoblaciones normales independientes, esta variable se distribuye en cada subpoblación con media y varianza diferentes. El resultado de esto es una mixtura dada.

El objetivo es estimar la media y la varianza de las subpoblaciónes, asi como la probabilidad de que un individuo pertenezca a alguna de las mismas.

Problema módelo

Es importante notar que no se tiene información de a que subpoblación pertenece cada individuo.

Población

Módelo

La densidad del módelo es:

\[ f(x;\mu_1,\mu_2,\sigma^2_1,\sigma^2_2,\pi)=\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2) \]

Con esto la función de log-verosimilitud queda:

\[ l(x;\mu_1,\mu_2,\sigma^2_1,\sigma^2_2,\pi)= \]

\[\sum_{i=1}^n \ln(\pi d(x_i;\mu_1,\sigma^2_1)+(1-\pi)d(x_i;\mu_2,\sigma^2_2)) \]

Ecuaciónes normales

Obteniendo un sistema a partir de las derivadas tenemos:

\[ \frac{dl}{d\pi}=\sum_{i=1}^n \frac{d(x_i;\mu_1,\sigma^2_1)-d(x_i;\mu_2,\sigma^2_2)}{\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2)} \] \[ \frac{dl}{d\mu_k}=\sum_{i=1}^n \frac{(-1)^{k+1}((3-2k)\pi+k-1)\frac{d}{d\mu_k}d(x_i;\mu_k,\sigma^2_k)}{\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2)} \]

\[ \frac{dl}{d\sigma^2_k}=\sum_{i=1}^n \frac{(-1)^{k+1}((3-2k)\pi+k-1) \frac{d}{d\sigma^2_k}d(x_i;\mu_k,\sigma^2_k)}{\pi d(x;\mu_1,\sigma^2_1)+(1-\pi)d(x;\mu_2,\sigma^2_2)} \]

El método de Newton

Sea \(f(x)\) una función de la cúal deseamos encontrar una raiz, sea \(\alpha\) esta raiz. Si realizamos la expansión de Taylor en el punto \(x_i\)

\[ f(x)=f(x_i)+f'(x_i)(x-x_i)+O((x-x_i)^2) \]

Si evaluamos en \(\alpha\) e ignoramos el segundo termino,

\[ 0=f(\alpha)=f(x_i)+f'(x_i)(\alpha-x_i) \] Así \[ \alpha=x_i-\frac{f(x_i)}{f'(x_i)} \]

Método de Newton

Este \(\alpha\) que estamos proponiendo no es exactamente la raiz (ya que aproximamos al descartar los terminos de orden 2), pero esperamos que este más cerca de la raiz, el método entonces esta caracterizado por la sucesión:

\[ x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)} \]

Esto se generaliza a funciones \(f:\mathbb{R^n}\to\mathbb{R^m}\) (Método de Newton-Rhapson) con

\[ \vec{x_{i+1}}=\vec{x_i}-J(\vec{x_i})^{-1}f(\vec{x_i}) \]

Método de Newton: Ventajas

  • Es aplicable para cualquier función (difenciable con algunas restricciones en la segunda derivada)
  • Converge de manera cuadratica cuando se cumplen supuestos de suavidad de la función (Aproximadamente duplica los digitos exactos en cada iteración)

Método de Newton: Desventajas

  • Puede no converger (converge en una vecindad de la solución).
  • Hay dificultad operacional al calcular la derivada (la matriz \(J\) y su inversa)
  • Encuentra solo una solución (en caso de existir varias)
  • No aprovecha la estructura del problema (optimización)

Paso 1: Dar un valor inicial a x.

Paso 2: Calcular la recta tangente de la función en el punto dado.

Paso 3: Usar la recta tangente para actualizar el x.

Paso 4: Repetir hasta la convergencia.

Visualización

Divergencia en el método de Newton

Ejemplo: Programación en R

func=function(x){(x^3-6*x^2+11*x-6)*1}
dfunc=function(x){(3*x^2-12*x+11)*1}
valor=list(-8)
for(i in 2:10){
        valor[[i]]=valor[[i-1]]-func(valor[[i-1]])/dfunc(valor[[i-1]])
}
data.frame(unlist(valor))
##    unlist.valor.
## 1     -8.0000000
## 2     -4.6889632
## 3     -2.4927804
## 4     -1.0454795
## 5     -0.1060077
## 6      0.4819019
## 7      0.8168003
## 8      0.9646914
## 9      0.9982722
## 10     0.9999955

Método EM (Expectation Maximization)

En el ejemplo inicial tenemos una variable que no fue observada, de la cual solo podemos conjeturar su valor.

  • Es posible estimar la esperanza de esta variable si tuvieramos los valores de los parametros de media y varianza de las subpoblaciones.
  • Es posible estimar los parametros de las dos distribuciones normales si tuvieramos la probabilidad de pertenencia a cada grupo.

Pasos del método

  1. Inicializar el proceso con algun valor para los valores de los parametros a estimar.
  2. Optimizar la verosimilitud de los parametros ocultos considerando que los parametros visibles son conocidos (este es el paso E, ya que es equivalente a tomar la esperanza bajo los parametros desconocidos).
  3. Optimizar la verosimilitud de los parametros visibles considerando que los parametros ocultos son conocidos (Este es el paso M).
  4. Repetir hasta la convergencia

Representacion Gráfica

Drawing

Expectation Maximization: Ventajas

  • Siempre converge.
  • Relativamente intuitivo y facil de implementar.

Expectation Maximization: Desventajas

  • Solo funciona para problemas de optimizacion que cumplen cierta estructura.
  • La velocidad de convergencia es baja.
  • Esta convergencia no se puede asegurar que sea a un optimo global.
  • Encuentra solo una solución (en caso de existir varias).

Para mixtura

En el caso especial de una mixtura (ejemplo inicial).

  • Inicie los parametros \(\mu_1,\sigma_1^2,\mu_2,\sigma_2^2.\pi\)
  • Expectation Step: calcule las probabilidates

\[ \gamma_i=\frac{\pi \Phi_{\mu_2}(y_i)}{\pi \Phi_{\mu_2}(y_i)+\pi \Phi_{\mu_1}(y_i)} \]

\[ \pi=\bar \gamma_i \]

Para mixtura

  • Paso de Maximización: Obtenga el estimador maximo verosimil para \(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2\). Que son los estimadores usuales tomando como pesos las probabilidades \(\gamma_i\) calculadas en el paso anterior.

  • Repita los pasos dos y tres hasta la convergencia

Paso 0.

Se utilizaron números aleatorios entre 0 y 10 para inicializar las medias, las varianzas se iniciaron cada una con la mitad de la varianza muestral.

Paso 1.

Se actualiza la probabilidad de que cada punto pertenezca a cada poblacion (es decir se calcula \(\gamma_i\))

Paso 2.

Con esta información se estiman de nuevo medias y varianzas usando el estimador maximo verosimil.

Paso 3.

Repetir hasta la convergencia

Visualización

Ejemplo: Programación en R

data(iris)
qplot(data=iris,x=Petal.Length,y=Sepal.Length,color=Species)

Ejemplo: Programación en R

library(mclust)
mod=Mclust(iris[,1:4])
summary(mod)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm 
## ----------------------------------------------------
## 
## Mclust VEV (ellipsoidal, equal shape) model with 2 components:
## 
##  log.likelihood   n df       BIC       ICL
##        -215.726 150 26 -561.7285 -561.7289
## 
## Clustering table:
##   1   2 
##  50 100

Ejemplo: Programación en R

pd=predict(mod)
qplot(data=iris,x=Petal.Length,y=Sepal.Length,
      color=as.factor(pd$classification))

Ejemplo: Programación en R

library(mclust)
mod=Mclust(iris[,1:4],G=3)
summary(mod)
## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm 
## ----------------------------------------------------
## 
## Mclust VEV (ellipsoidal, equal shape) model with 3 components:
## 
##  log.likelihood   n df       BIC       ICL
##       -186.0736 150 38 -562.5514 -566.4577
## 
## Clustering table:
##  1  2  3 
## 50 45 55

Ejemplo: Programación en R

pd=predict(mod)
qplot(data=iris,x=Petal.Length,y=Sepal.Length,
      color=as.factor(pd$classification))

Ejemplo: Programación en R

table(as.factor(pd$classification),iris$Species)
##    
##     setosa versicolor virginica
##   1     50          0         0
##   2      0         45         0
##   3      0          5        50

Referencias

Gracias

  • Esta presentación se realizó usando RMarkdown (ioslides)
  • Las gráficas de esta presentación se realizaron usando ggplot2
  • Las animaciones se realizaron usando plotly
  • Esta presentación y el código fuente están en https://github.com/kuhsibiris/Newton-y-EM bajo licencia Apache v2.0